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MARCHING  ITERATIVE  METHODS 

FOR  THE  PARABOLIZED  AND  THIN  LAYER  NAVIER-STOKES  EQUATIONS 


Moshe  Israeli 

Technion  -  Israel  Institute  of  Technology,  Haifa  (Israel) 


Abstract 


Downstream  marching  iterative  schemes  for  the  solution  of  the 
Parabolized  or  Thin  Layer  (PNS  or  TL)  Navier-Stokes  equations  are 
described.  Modifications  of  the  primitive  equation  global  relaxation 
sweep  procedure  result  in  efficient  second-order  marching  schemes. 

These  schemes  take  full  account  of  the  reduced  order  of  the  approximate 
equations  as  they  behave  like  the  SLOR  for  a  single  elliptic  equation. 
The  improved  smoothing  properties  permit  the  introduction  of  Multi-Grid 
acceleration.  The  proposed  algorithm  is  essentially  Reynolds  number 
independent  and  therefore  can  be  applied  to  the  solution  of  the  subsonic 
Euler  equations.  The  convergence  rates  are  similar  to  those  obtained  by 
the  Multi-Grid  solution  of  a  single  elliptic  equation;  the  storage  is 
also  comparable  as  only  the  pressure  has  to  be  stored  on  all  levels. 
Extensions  to  three-dimensional  and  compressible  subsonic  flows  are 
discussed.  Numerical  results  are  presented.  .  1  J  ‘ 
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1.  INTRODUCTION 


Considerable  evidence  accumulated  recently  about  the  applicability 
of  the  Parabolized  Navier-Stokes  equations  for  high  Reynolds  number  flows 
with  a  principal  flow  direction;  see  Rubin  (1).  The  PNS  equations  are 
obtained  by  neglecting  the  streamwise  viscous  terms  in  the  Navier-Stokes 
(NS)  equations.  When  the  viscous  terms  in  the  circumferential  direction 
are  also  neglected,  one  gets  the  Thin  Layer  approximation. 

The  steady  PNS  equations  still  have  an  elliptic  nature,  and  there¬ 
fore  the  initial  value  problem  in  the  marching  direction  is  not  well 
posed  [2] .  A  well  posed  initial-boundary  value  problem  can  be  formulat¬ 
ed  by  specifying  (for  example)  upstream  and  side  conditions  for  the 
velocities  and  one  downstream  condition  for  the  pressure.  This  coupled 
system  of  partial  differential  equations  behaves  like  a  single  elliptic 
equation  for  the  pressure.  Therefore  the  PNS  equations  must  be  solved 
globally  and  cannot  be  solved  by  a  single  sweep  marching.  The  reduced 
order  of  the  PNS  equation  can  be  exploited  by  constructing  an  iterative 
marching  method  for  updating  the  pressure  field  only.  Such  a  multiple 
sweep  iteration  method  has  the  advantage  that  the  velocity  fields  are 
generated  during  the  marching  process  and  only  the  pressure  field  has  to 
be  stored  from  sweep  to  sweep.  A  considerable  saving  in  storage  results. 
However,  simple  minded  marching  does  not  result  in  good  convergence 
properties  and  sometimes  diverges.  For  the  two-dimensional  incompres¬ 
sible  case,  Israeli  and  Lin  [3]  devised  a  stable  marching  scheme  that 
behaves  like  the  Successive  Line  Over  Relaxation  (SLOR)  method  for  a 
single  elliptic  equation.  The  good  smoothing  properties  of  the  above 
mentioned  scheme  can  be  used  in  a  Multi-Grid  (MG)  framework  in  order  to 
accelerate  the  convergence  of  the  solution  of  the  PNS  (or  TL)  equations. 
The  marching  scheme  is  implemented  using  a  new  stable  algorithm  which  is 
second  order  also  in  the  marching  direction.  The  same  method  can  be 
used  without  modification  for  the  subsonic  Euler  equations  as  the  effect 
of  the  Reynolds  number  on  the  convergence  rate  is  insignificant.  In  two 
dimensions  the  PNS  and  TL  equations  are  identical  and  therefore  the  same 
analysis  applies  to  both. 
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It  turns  out  that  the  extension  to  three-dimensions  is  conceptually 

simple;  but  the  resulting  algorithm,  a  successive  plane  over  relaxation, 
is  complicated  by  the  requirement  of  the  simultaneous  solution  of  the 
equations  in  planes  perpendicular  to  the  marching  direction.  This 
problem  can  be  alleviated  by  splitting  of  the  equation  of  continuity 
from  the  momentum  equations. 

The  extension  of  the  method  to  compressible  flows  is  conceptually 
non  trivial.  The  original  iterative  method  is  based  on  the  concept 
that  the  convergence  relies  on  the  implicit  relaxation  of  a  single 
quantity,  the  pressure,  which  approximately  satisfies  a  single  elliptic 
equation.  In  the  compressible  case  a  viable  approach  is  to  eliminate 
the  pressure  and  to  derive  an  equation  for  p,  the  logarithm  of  the 
density.  It  can  be  shown  that  p  satisfies  approximately  equation  (1.1), 

(1  ‘  m2)P98  +  Pnn  “  °*  (K1) 

where  M  is  the  Mach  number  and  s  and  n  are  coordinates  along  and 
perpendicular  to  the  flow  direction.  Although  this  equation  la  never 
derived  or  used  In  the  algorithm,  It  reveals  the  fact  that  for  M  <  1 
the  upstream  influence  is  transmitted  through  the  quantity  p,  and 
therefore  only  this  quantity  should  be  stored  or  updated.  The  flow  of 
Information  should  be  downstream  for  the  velocity  and  temperature  and 
upstream  for  the  density,  and  the  difference  scheme  must  be  built 
accordingly.  For  supersonic  flows  the  flow  of  information  should  be 
only  downstream  and  the  marching  method  is  non-iterative.  For  super¬ 
sonic  flows  with  Imbedded  subsonic  regions,  the  iterative  method  should 
be  used,  combined  with  an  appropriate  switching  at  shock  waves  and  sonic 
lines . 

It  should  be  pointed  out  here  that  the  present  approach  is  very 
different  conceptually  from  that  of  Reddy  and  Rubin  (4).  Although  they 
used  our  idea  (Israeli  (6])  of  back  shifting  the  pressure,  one  full  mesh 
distance,  with  respect  to  the  velocity  for  incompressible  flows,  their 
generalization  to  compressible  flows  is  a  Mach  number  dependent  shift 
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which  vanishes  for  M  >  1.  This  smooth  transition  from  subsonic  to 
supersonic  flows  is  questionable  since  the  change  of  type  of  equation 
(1.1)  is  sudden,  at  M  •  1.  Indeed,  only  our  full  6hift  is  used  in  their 
papers  and  properly  results  in  a  conservative  scheme  across  a  shock. 

Another  question  raised  by  the  above  mentioned  paper  is  that  of  the 
distinction  between  the  pressure  which  uses  downstream  data  and  the 
density  which  uses  upstream  data.  This  obscures  the  issue  of  the  direc¬ 
tion  of  flow  of  information  and  proper  location  of  boundary  conditions. 
This  approach  should  result  in  inconsistency  of  boundary  data  and  may 
eventually  lead  to  ill  posedness  and  divergence. 

In  the  next  sections  we  will  summarize  our  previous  theoretical 
results,  present  some  new  numerical  results  and  the  extensions  to  3-D 
and  compressible  flows. 


2.  FORMULATION  FOR  THE  INCOMPRESSIBLE  CASE 

For  simplicity  we  will  consider  initially  the  case  of  the  steady, 
incompressible,  and  two-dimensional  PNS  (or  TL)  equations  in  cartesian 
coordinates  lx,-y]  : 


U  +  V  =  0 
x  y 


(U  )  +  (UV)  «  -P  +  U  /Re 

x  y  x  yy 


(UV)  +  (V  )  =  -P  +  V  /Re 

x  y  y  yy 


(2.1) 


(2.2) 


(2.3) 


where  x  is  the  mainstream  direction.  Re  is  the  Reynolds  number.  U 
and  V  are  the  nondimensional  velocity  components  in  the  x  and  y 
direction,  respectively.  P  is  the  nondimensional  pressure. 


The  two-dimensional  NS  equations  are  elliptic  of  order  four  -  Brandt 
and  Dinar  IS] .  The  PNS  are  elliptic  only  of  order  two  like  the  Poisson 
equation  (the  mathematical  nature  of  several  two-dimensional  and  three- 
dimensional  approximations  to  the  Navier-Stokes  equations  was  analyzed 
in  17]).  This  ellipticity  is  due  to  the  pressure  gradient  terms  via 
the  continuity  equation.  A  well  posed  problem  can  be  formulated  by 
defining  the  boundary  conditions  as  described  in  Fig.  2.  The  following 
Dinchlet  conditions  may  be  specified: 


v’ls'.v 


*  upstream  boundary  (AB) s 

*  at  a  solid  wall  (AD) : 

*  at  the  outer  boundary  (BC) : 


Other  boundary  conditions  can  be  used*  but  the  same  number  of  conditions 
on  each  boundary  must  be  kept. 

In  order  to  separate  linear  and  non  linear  effects,  some  of  the 
convergence  tests  were  performed  with  the  following  linear  version  of 
equations  (2.1)- (2.3) : 


1 .  ? 
in 

V  « 

Vin 

(2.4) 

wall  * 

V  «= 

Vwall 

(2.5) 

out 

V  * 

Vout 

(2.6) 

* 

down* 

(2.7) 

U  +  V 
x  y 


(aU)  +  (bU)  -  -P  +  U  /Re 
x  y  x  yy 


UV)  +  (bv)  =  -P  +  V  /Re 

x  y  y  yy 


(2.8) 

(2.9) 

(2.10) 


where  a  and  b  are  known  functions  of  x  and  y. 

3.  DISCRETIZATION  AND  MARCHING 

Numerical  solutions  of  Eqs.  (2.1)— (2.3)  are  obtained  by  spreading 
a  grid  over  the  computational  domain.  Let  us  assume  that  the  grid  points 
are  distributed  evenly  along  the  x  and  y  coordinates  with  the  spac¬ 
ing  Ax  and  Ay  respectively.  When  differencing  these  equations  it 
should  be  remembered  that  their  nature  should  be  reflected  (1,B)  in 
the  finite  difference  approximation.  In  order  to  be  consistent  with  the 
boundary  layer  (parabolic)  nature  of  the  flow,  the  axial  gradients  of 
the  velocities  should  be  computed  using  only  upstream  values,  while  the 
elliptic  nature  is  preserved  by  forward  differencing  the  axial  pressure 
gradient  11,8,9).  Consequently,  it  was  assumed  that  a  stable  marching 

scheme  must  be  of  the  first  order  in  the  marching  direction.  It  turns 
out  that  this  effect  can  be  achieved  by  a  judicious  choice  of  the  place¬ 
ment  of  the  variables  to  be  solved  at  each  station.  The  choice  can  be 
explained  most  easily  by  taking  V  =  0  and  ^-  =  0  in  Eq.  (2.2)  for  U, 
yielding 


A.V.Vj.V.  aVVa’  ■  ■  o  • 


U  *  -  p  . 

X  X 

A  first  order  difference  scheme  then  becomes 
2  2 

m,}  m-l,j  m,j  *m+l,g  • 

the  unknowns  are  U  and  p  . .  The  scheme  first  suggested  by  Israeli 
m,3  *m,;j  ”  1 

[9,10]  is: 

u2  -  u2  =  p  -  p 

m+l,j  m,j  m,j  m+l,j 

with  the  unknowns  U  ,  .  and  p  ..  The  scheme  is  centered  about  m+  — 

m+l,]  m,3  2 

and  is  second  order.  This  approach  was  subsequently  used  by  Rubin  and 
Reddy  [8]  and  Reddy  and  Rubin  [4J. 


In  addition,  one  may  stagger  the  velocity  V  with  respect  to  the 
other  variables  as  shown  in  Fig.  3,  where  the  centering  points  of  the 
different  difference  equations  are  also  plotted.  The  differential  equa¬ 
tions  are  approximated  by  central  second-order  approximations  whenever 
needed  averaging  was  used  as  is  usually  done  for  staggered  grids. 

Numerical  experiments  with  a  first  order  computer  code  show  that 
the  solution  after  one  marching  sweep  is  not  close  to  the  final  solution 
of  the  PNS  equations  when  the  initial  pressure  field  is  constructed 
using  the  boundary  layer  assumption  p  “0.  Since  the  p  term  is 

y  * 

forward  differenced,  some  global  iterations  over  the  whole  solution 
domain  should  be  performed  in  order  to  converge  the  explicit  contribu¬ 
tion  to  this  pressure  terjn.  The  simplest  global  iterative  technique  to 
solve  the  equations  is  by  multiple  marching  sweeps  with  the  primitive 
equations  where  only  the  pressure  field  is  kept  from  iteration  to  itera¬ 
tion  [1].  Numerical  experiments  also  show  that  for  certain  nets  this 
procedure  diverges.  The  divergence  occurs  also  for  the  linearized  ver¬ 
sion  of  Eqs.  (2.1)  -(2.3).  Figure  1  presents  the  residual  of  the  pressure 
field  as  a  function  of  the  global  iteration's  sweep  number  for  a  21  *  11 
field.  A  jump  is  encountered  every  10  iterations  (probably  related  to 
the  arrival  of  the  boundary  pressure  pulse  traveling  at  the  numerical 
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scheme  speed)  leading  to  ultimate  divergence.  However  convergence  was 
reported  with  different  mesh  and  boundary  conditions  and  also  when 
combining  the  above  procedure  with  a  multigrid  technique  [4],  It  was 
thought  that  the  replacement  of  one  of  the  momentum  equations  by  the 
Poisson  equation  for  the  pressure  will  improve  the  convergence  rate,  but 
the  solution  did  not  satisfy  the  replaced  momentum  equation.  A  success¬ 
ful  implementation  of  the  marching  technique  is  derived  in  the  next 
section.  A  short  and  reduced  version  of  the  analysis  was  presented  first 
in  J3). 


4.  A  MULTI -GRID  ALGORITHM 

The  Multi-Grid  technique  is  a  numerical  strategy  for  substantially 
improving  the  convergence  rate  of  an  iterative  procedure.  In  order  to 

facilitate  comparison  with  theory  ,  the  accomodative  C-cycle  MG  algorithm 

* 

was  chosen. 

Each  MG  process  consists  of  three  basic  parts:  relaxation,  restric¬ 
tion,  and  interpolation  [5]. 

The  Relaxation  Scheme 

The  overall  convergence  rate  of  any  MG  process  is  greatly  influenc¬ 
ed  by  the  smoothing  properties  of  the  relaxation  scheme.  It  can  be  shown 
analytically  and  experimentally  that  the  usual  multiple  sweep  marching 
[1]  does  not  have  good  convergence  and  smoothing  properties  because  short 
wave  errors  are  not  efficiently  smoothed.  Israeli  and  Lin  [3]  showed 
that  certain  modif ications  in  the  streamwise  momentum  equation,  which 
vanish  upon  convergence,  give  rise  to  an  iterative  scheme  which  is  equi¬ 
valent,  in  the  linear  case,  to  the  SLOR  method  for  one  Poisson  equation. 

In  the  general  nonlinear  case  the  modified  iterative  process  is  essential¬ 
ly  equivalent  to  the  relaxation  of  a  single  nonlinear  Poisson-like  equa¬ 
tion  for  the  pressure.  The  velocities  can  be  viewed  as  auxiliary  vari¬ 
ables  needed  during  the  marching  since  they  have  no  "memory"  by  them¬ 
selves. 

Fur thermore ,  we  have  automatically  gained  the  good  smoothing 
properties  of  the  line  relaxation  scheme  of  a  single  Poisson  equation. 

The  problems  associated  with  the  loss  of  ellipticity  of  the  difference 

*  Some  of  the  elements  of  the  present  approach  were  used  independently  by 
Rubin  and  Reddy  18).  Detailed  comparisons  cannot  be  made  because  converg¬ 
ence  rates  and  storage  estimates  were  not  presented  there. 


approximation  for  the  Navier-Stoke?  equations  at  high  Reynolds  number 
( 5 J  are  thus  avoided  and  no  upstream-weighting  or  artificial  viscosity 
are  required.  There  results  a  considerable  saving  in  storage,  as  well 
as  a  simpler  relaxation  scheme  (compare  to  the  distributive  relaxation 
[5])  where  the  convergence  rate  is  essentially  independent  of  the 
Reynolds  number.  We  note  that  the  same  marching  algorithm  can  thus  be 
used  for  the  (subsonic)  Euler  equation  with  the  same  favorable  converg¬ 
ence  rate.  (For  supersonic  flows  the  marching  method  is  non  iterative.) 

A  part  of  the  analysis  of  (3)  is  repeated  here  to  motivate  the 
later  extensions  to  three-dimensional  and  compressible  flows.  We  start 
with  the  PNS  equations  (2.1) -(2. 3)  and  linearize  them  about  a  constant 
state.  We  also  introduce 

2 

L(f)  =  q— —  (Uf)  +  (Vf)  -  2__  f  (4.1) 

9x  3y  Re  _  2 

3y 

where  U  and  V  are  constant  reference  velocities.  The  next  step  is 
to  discretize  the  equations  only  in  the  x  direction  to  obtain: 


U 

m-1  m 


Ax 


+  (V  )  =0 

y  m 


(4.2) 


P  -P 

D(U  ) - Eli-E 

m  Ax 


(4.3) 


D  (V  )  =  -  (P  ) 
m  y  m 


(4.4) 


where  U  .  V  and  P  are  functions  of  y.  Here  D(f  )  is  the  semi 
m  m  m  m 


discretized  form  of  L(f)  at  the  marching  station  m.  The  semi  - 
discretized  system  should  be  discretized  also  in  the  y  direction  before 
solution  is  attempted,  but  since  the  specific  form  of  this  discretization 
is  not  important  for  the  following  argument,  we  postpone  this  step  for 
the  sake  of  transparency. 


k  k 

The  marching  iterative  procedure  assumes  that  U  (y) ,  V  (y)  are 
k- 1 

known  as  well  as  P  for  m  =  2,3,4,...  M,  where  k  is  the  current 

m 


iteration  index.  Therefore,  the  mai  china  r.-i.eme  for  m  5  2  is: 


k  k 
U  -U 
m-1  m 


+  (V  ) 


0 


(4.5) 


\  \  V  v\  *.  \  •-  ■.  v  ■ ».  «  '  *_ 
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D(Vk)  =  -  (Pk)  . 
m  y  n 


We  now  apply  D  to  Eq.  (4.5)  and  differentiate  Eq.  (4.7)  with 
respect  to  y.  Elimination  of  the  V  terms  between  Eqs.  (4.5)  and 
(4.7)  gives: 


D(Uk  ,  -  Uk)  =-(Pk  )  Ax  . 
m-1  m  yy  m 

Now  by  substitution  of  Eq.  (4.6)  into  Eq.  (4.8)  we  get: 


k-1  k  k-1  k  .  2,  k  , 

P  -  P  -  P  +  P  ,  +  Ax  (P  )  =  0 

m+1  in  m  n-1  yy  m 


(4.9) 


It  follows  that  the  marching  scheme  for  the  primitive  system  (4.2)- (4.4) 
can  be  viewed  as  a  line  iterative  scheme  for  the  semi-discretized  Laplace 
equation;  indeed  upon  convergence  Eq.  (4.9)  will  become: 


P  -2P  +P  , 
m+1  m  m-1 


+  (P  )  =0 

yy  m 


(4.10) 


In  order  to  find  out  the  rate  of  convergence  of  Eq.  (4.9)  to  the  final 
state  (4.10),  we  Fourier  transform  Eq.  (4.9)  in  y  assuming  appropriate 
boundary  conditions  in  that  direction: 


k  -  I  ny  k-1  _  I  ny 

P  =  z  e  J  ,  P  =Ze 
n  m  mm 


(4.11) 


where  I  =  -1  and  n  is  the  Fourier  wave  number.  After  substituting 
these  definitions  into  Eq.  (4.9)  we  get: 


2  2- 

Z  -  Z  -Z  +Z  -AxnZ  =0. 
m+1  m  m  m-i  m 


(4.12) 


Transforming  in  the  x  direction  we  define  again 


Z  =  Ae 


i  ex 


(4.13) 


where  0  is  the  wave  number  in  the  x  direction.  By  substitution  of 
this  definition  into  Eq.  (4.12)  we  get: 


(4.14) 


1  ,  .22  -10 
1+Ax  n  -e 


This  means  that  all  the  long  waves  (with  small  Ax  n  )  in  the  cross  flow 


direction  are  only  weakly  damped  irrespective  of  their  structure  in  the 
marching  direction.  In  particular  the  n=0  modes  which  exist  for 
derivative  boundary  conditions  in  the  cross  flow  direction  are  not 
affected  at  all  by  the  relaxation,  i.e.. 


On  the  other  hand,  the  well  known  SLR  scheme  for  the  Poisson  equation 
gives  (after  the  same  Fourier  transformations) : 


Z  -  (2+  Ax2n2)Z  +  Z  =0 
m+1  m  m-1 


and 


I- 

'a 


q-e 


-16 


2  2 

•  q  =  2  +  Ax  n  )  2 


and  also 


i  A ,  2 


'A1  2 

q  +l-2qcos8 


(4.15) 


(4.16) 


(4.17) 


This  quantity  is  less  than  1  for  all  acceptable  q's  and  cos0  <  1. 

Most  waves  are  strongly  damped,  and  only  the  longest  waves  in  both  direc 
tions  are  weakly  damped  by  the  iteration.  This  behavior  was  used  to 
accelerate  the  convergence  as  is  done  by  the  SLOR  technique,  Chebychev 
acceleration,  or  Multi-Grid  method.* 

The  question  is  how  to  generate  an  equivalent  relaxation  scheme 
for  the  primitive  system  in  the  marching  form.  This  means  that  we  may 
add  terms  which  can  be  evaluated  during  the  marching  process  but  should 
vanish  upon  convergence. 

A  rational  approach  to  the  construction  of  the  relaxation  scheme 
is  to  retrace  backwards  the  steps  of  the  derivation  of  the  discrete 
Laplace  equation  from  the  discrete  primitive  equations.  We  start  from 
the  SLOR  equation  (4.15): 


Z  .  -  2Z  +  Z  . 
m+1  m  m-1 


.  2  2- 
Ax  n  Z 


which  we  inverse  Fourier  transform  with  respect  to  y  to  get: 


P*-1 

m+1 


-  2P  +  P 
m 


m-1 


-Ax2(Pk  ) 
yy  m 


*It  was  pointed  out  by  J.  South  that  (4,3)  can  be  viewed  as  an  over¬ 
relaxed  version  oi  (4.10)  with  an  over-relaxation  factor:  0=2  which 
is  not  a  good  *»hoice  for  0. 


)  / 


Nov.’,  we  substitute  Eq.  (4.8)  for  the  right  hand  side  of  the  last  equa¬ 
tion  to  get: 


_k-l  _k  k 

P  -  2P  +  P 
m+1  m  m-1 


which  can  be  written  as: 


_k  ,k-l 

P  )  +  (P 
m  m 


ixD(Uk  ,  -  Uk) 
m-1  m 


PK) -  (Pk_1 -  Pk  , )  =  AxD (u  -  u  )  ; 

m  m  m-1  n-1  m 


adding  the  equations  from  m  =  2  and  using  the  linear  form  of  D  we  get: 


.  k-1  k.  r  .  K-i  k,  ,  k-1  ^k, 

-(P  -  P  )  -  )  (P.  -P.)+(P.  -  P,  )  =  AxD  (U  -U, 

m+1  m  .  *•_  a.  l  2  1  ml 


k.  ,  k-1  V, 


but  from  Eq.  (4.3) 


P2~1  "  P1  =  ”AxDtu^); 


(4.18) 


therefore,  we  get  for  m  >  2 


k-1  k 

D(u  ) - Eii-E  -  i-  I  IP11-1 

m  Ax  Ax  .*■_ 

1=2 


(4.19) 


Eq.  (4.19)  contains  all  the  modifications  required  in  order  to 
convert  the  iteration  scheme  of  (4.5)— (4.7)  into  a  scheme  equivalent  to 
the  SLOR  scheme  for  one  Laplace  equation  with  "over-relaxation"  factor 
(1=1.  We  see  that  in  this  approach  only  the  x  momentum  equation  is 
modified.  The  new  added  terra  can  be  generated  easily  during  the 
marching  process  and  is  inexpensive  in  storage  (one  extra  line  vector) 
and  computation  (one  substruction  per  grid  point).  In  what  follows  we 
will  derive  Eq.  (4.19)  in  a  more  general  way  and  introduce  the  over- 
relaxation  parameter  (1  >  1. 

In  practice  we  will  use  difference  approximations  and  boundary 
conditions  also  in  the  y  direction,  and  the  resulting  scheme  may  not  be 
amenable  to  the  discrete  analogue  of  the  Fourier  transform.  It  is 
therefore  worthwhile  to  generalize  the  previous  approach  by  using  the 
matrix  finite  difference  formulation. 

Let  the  vectors  U  ,  V  ,  P  contain  the  N  values  of  the  cor- 
rn  m  rn 

responding  variables  on  the  m-th  line  (x  =  constant)  of  the  marching 
sweep  (including  the  specified  boundary  values).  The  p-momentum 


,\  .'w  >  v  v  ■■  ■  ■  *  .•  *  i  ■  j  v.  ■ .  *  .*,■ 


equation  (1.5)  can  be  written  in  the  form: 


P  -  P  =  -AxD (U  )  =  R  . 

m+1  m  m  m 


(4.20) 


On  the  other  hand  elimination  of  V  between  the  continuity  and  the 

m 

V-momentum  equations  will  result  in: 


FP  =  R  -  R  . 
m  m  m-1 


(4.21) 


where  F  =  Ax  I  — y  .  Substructing  successively  u-momentum  equations 
(4.20)  and  using  Eq.  (4.21)  gives: 


P  ,  -  HI+F1P  +  P  ,  =  0,  m  =  3,4,.. 
m+1  m  m-1 


(4.22) 


which  is  Laplace's  equation.  The  first  equation  of  (4.20)  can  be  used 
as  a  derivative  condition  at  the  left  (inlet)  boundary,  namely: 


P3  -  P2  =  R2 


(4.23) 


We  now  apply  the  SLOR  scheme  to  the  last  two  equations  (ignoring 
temporarily  the  downstream  boundary  condition)  to  get  the  downstream 
marching  form: 


+  p(k-l>  =  R, 

2  3  2 


P(k*  -  (21  +  F )  P  +  =0,  m  =  3,4,.. 

m-1  m  m+1 


(4.24) 


(4.25) 


(k  1  *  (k-1 ) 

where  P  =  OP  +  (1-0)P  :  f 1  is  the  overrelaxation  factor,  and 

mm  m 

the  superscript  denotes  the  iteration  sweep  number.  In  order  to  recover 
the  primitive  variable  formulation,  we  relate  the  velocity  field  in  Eq. 
(4.21)  to  the  starred  pressure  field,  i.e., 


FP  *  R  -  R  .  . 
m  m  m-1 


(4.26) 


Substitution  in  Eq.  (4.25)  gives: 


P(k)  -  2P*  +  p(k-1)  =  r  -  R  ,  m  =  3,4,... 
m-1  m  m+1  m  m-1 


Successive  summations  of  Eqs.  (4.24)  and  (4.25)  giv6. 


(4.27) 


(k-1) 


P  =  R  +  S  , 
m  m  m 


m  =  2,3,4,... 


(4.28) 


equation . 


which  1 £  the  primitive  variable  marching  form  of  the  u-momentum 

The  source  term:  S  in  £q .  (4.26)  satisfies: 

m 


*  (k-1 )  * 

S  +  (P  -  PU  )  +  (F- 
m-i  m  m  m-1 


r  '  )  ,  m  =  3,4,... 
m- 1 


(4.29) 


with  S  =  0.  It  can  be  seen  that  S  vanishes  upon  convergence  The 
2  m  r 

computational  form  of  (4.28)  for  m  =  3,4,...  is: 


-2P  =  R  +  S 

m  m  m 


(4.30) 


s  =  5  -  p(k-1}  4  2P*  -  P(k>-  c  -  _p*  .  p* 

m  m-1  m4l  m-1  m-1'  "2  P2  P3 


(4.31) 


Thus,  the  theory  of  overrelaxation  can  be  applied  exactly  to  the 
constant  coefficient  case  of  system  (2 . 8) -  (2 . 10) .  For  the  non-linear 
case  this  theory  can  serve  as  a  guide  to  the  choice  of  (2.  Alternately, 
one  can  choose  (2=1  and  apply  the  Multi-Grid  procedure. 

Restriction  and  Storage  Requirements 

Let  the  finite  difference  approximation  of  equations  (2.1)- (2.3) 
on  the  finest  grid  M  be  represented  as  in  [5] : 


LMWM(x)  =  FM(x) 


(4.32) 


-  - M  M  M  M  T 

where  x  =  (x,y) ,  w  =  [U  ,V  ,P  ]  is  the  exact  solution  of  the  dif¬ 
ference  equations,  and  j  is  the  number  of  the  differential  equation, 
j  =  1,2,3. 

The  problem  is  transferred  from  the  current  level  k  to  a  coarser 
level  k-1,  see  Fig.  4,  by  correcting  the  right  hand  side  of  (4.32) 


_k-l,~,  . k-1  -k-l'k , ,  k-l,„k,~,  k-k,~.. 

F  .  (x)  =  L  .  (I  .  ,  v  (x)  )  *  I  .  ,  IF  .  (x)  -  L  .w  (x)  ] 

3  3  3  »k  3 ,k  3  3 


(4.33) 


in  the  Full  Approximation  Storage  (FAS)  mode,  w  (x)  is  an  approximation 

~k  -  k-1  -k-1 

to  Vf‘(x)  in  the  finer  level.  I.  ,  and  I.  ,  are  proper  restriction 

3 /k  3  ,k 

operators  for  equation  j. 

The  term  in  square  bracket  in  equation  (4.33)  is  the  residual  of  the 
j-th  equation.  For  the  present  marching  scheme  there  is  no  residual  in 
the  continuity  and  in  the  y-momentum  equations  since  they  are  solved 
exactly  in  each  step.  The  residual  of  the  x-momentum  equation  results 


only  from  the  Ftreamwise  pressure  Gradient  tern,  and  its  computation  need 

'k-1 

only  one  r abstraction .  I.  was  chosen  to  be  linear  interpolation, 

3  '  _ j  ~K— l-l-i  v - 1 

which  yields  for  the  continuity  ecuation:  L,'  (1,  w’lx))  =  0 .  I  , 

1  1  ,  V.  3  ,  k 

j  =  1,2  is  computed  by  averaging  in  both  the  x  and  y  directions. 

Ik  ^  is  a  simple  injection. 

In  summary,  equation  (4.33)  takes  the  following  terms: 

Fk_1(x)  =  0 

F*_1(x)  =  L^uli'Jw'tx))  +  Ik_1(Fk  -Lkwk(x)) 

2  2  2 ,  k  2 ,  k  2  2 

_k-l  .- .  k-1 . rk-1 -k  ' 

F  (x)  =  L  (I-,,w  ( x )  )  - 

3  3  3  ,k 

Two  consequences  should  be  emphasized: 

k-1  '  k-1  - 

(a)  Only  two  corrections  (F^  <x)  »  ^3  )  have  to  be  computed  and 

stored . 

(b)  All  the  dependent  variables  must  be  transferred  in  order  to  compute 

-k-l'k  - 

the  corrections  (L.  (I w  (x) )  ,  j  =  2 , 3)  .  Since  only  the  pressure 
3  D  <  k 

is  stored,  these  corrections  must  be  computed  during  the  marching 
process . 

It  follows  that  in  addition  to  the  pressure  on  all  grids,  one  has 
to  save  one  correction  term  for  each  momentum  equation  on  the  coarser 
grids.  Assuming  N  computational  points  on  the  finest  grid,  a  simple- 
minded  estimate  gives  35N/7  storage  locations  for  the  three-dimensional 
NS  Multi-Grid  solution,  and  11N/7  for  the  PNS  marching  MG  solution. 

For  the  two-dimensional  case  the  corresponding  figures  are  14N/3  and 
6N/3. 


Interpolation 


Since  the  present  marching  scheme  generates  the  velocity  field  from 
the  pressure,  only  the  correction  to  the  pressure  must  be  interpolated 
back  to  the  fine  grid. 


>.■  ■  '  -v" «-• \" ■v.'XTR'XT 


f. .  GENERALISATION’S 

Iri  ort-.cr  to  generalize  the  precedin';  approach  we  note  that  the 
essence  of  the  relaxation  procedure  is  tin  rrpl acomcnt  of  the  terr. 

Lx3P/3x  by  the  marching  difference  form: 

AxIr  =  -(2IC  +  V  (5-1} 

where  S  is  already  known.  If  (5.1)  is  differenced  it  will  (using 

the  definition  of  S  )  give  rise  to 
m 

*  *  -  -  t-i  '*  i- 

2  (P  -  P  ,)+S  -  S  =  -(P  .  -  2P  +  P‘  )  . 

m  m-1  m  m-1  m+1  m  m-1 

Thus,  the  correct  successive  line  over  relaxation  form  is  implicitly 
obtained  for  the  second  derivative  of  the  pressure  in  the  marching 
direction.  (It  should  be  emphasized  again  that  the  second  order  elliptic 
equation  for  the  pressure  is  neither  derived  nor  used  in  the  algorithm 
itself . ) 

The  implication  of  the  present  technique  is:  if  it  is  known  that 
the  equations  can  be  manipulated  so  that  some  variable  will  satisfy 
approximately  a  second  order  elliptic  equation,  we  should  use  the  replace¬ 
ment  (5.1)  for  the  derivative  of  that  variable  in  the  main  flow  direction. 
An  efficient  marching  scheme  will  thus  be  generated. 

The  present  version  of  the  algorithm  will  be  applied  to  the  sub¬ 
sonic  compressible  multi-dimensional  Navier-Stoke ’ s  equations.  Several 
particular  cases  will  be  examined. 


The  first  step  is  the  derivation  of  an  elliptic  equation  starting 

with : 

V-VV=--  Vp  +  vV2V  +  ^  V  ( v  »V)  (5.2) 

P  3 

where  V  is  the  velocity  vector.  In  addition  we  will  require  the  equa¬ 
tion  of  state  of  a  perfect  gas 


p  =  pP.T 


and  the  continuity  equation  in  the  form 


7»V  =  -V'V&np. 


It  follows  from  (5.3)  that 


lo  \i]  str  ear 


d  *■  p  ci  *■  n 

(1-M‘)  — *-  +  — *-  =  other  term. 

3s^  i).l‘ 

After  the  parabolization  of  the  viscous  tent,  only  ♦he  left  hand  rid 
has  a  second  derivative  of  p  m  the  rtreamwise  direction. 
Specifically  v7*-V*V  is  replaced  by 

v  — —  7  •  V  =  -v  — —  v*7in..  . 
er?  'dr\z 

This  term  cannot  become  large  since  the  pressure  does  not  have  larq 
gradients  in  the  boundary  layer. 

We  argue  that  if  our  iteration  is  appropriate  for  the  p  equa 
it  will  be  a  good  scheme  overall. 

To  net  a  successive  lir.e  (or  plane)  over  relaxation  scheme,  all 

C*  l  ) 

have  to  do  is  replace  all  the  occurrences  of  — ^  with,  the  mar  chine 

ex 

(S.l).  All  the  properties  of  Section  4  w.ll  be-  the  same  as  lone:  as 
M:  <  1. 

In  fact,  better  convergence  car:  be  expected  as  M-  approaches 
since  the  quantity  q  of  (4.  if.)  will  become  now  2+  (h>r  rm /1-K*  )  . 
Only  p  will  have  memory  and  must  be  globally  saved  and  updated  by 
iteration  procedure.  p  will  also  transmit  tiie  downstream  inform, at: 
and  must  be  specified  there. 

Tor  transonic  flows  a  cons  'rvat  i  or.  form  is  preferred  and  it  may 
mere  convenient  to  work  with  r  rather  than  Inc.  A:,  clirptic  t 
car.  be  derived  for  c,  but  care  rust  be  taken  to  transmit  t  ne  uowr.st  r 
information  via  c  .  Upstream.  :  nformi  it  ion  srcul  d  net  be  ♦  rarer:  tt 
p  at:  2  p  siioulo  not  be  rr1,  ifjou  ,»t  the  inflow;  ot.uc-rwi « •  ,  tin  rr. 

for  t-xnm;  1  e  , 


v>u  1 1  be  over :  p.eci f  i ed  . 
discretized  as 


Coi.'  i  dor  , 


t  b  e  t. r '  r  rr 


I  K  vi  y  .  w  .1 

v  m  m  m- 1  m-1 ' 

however,  at  the  station  tn  we  should  compute  the  Um  velocities  coupled 
with  the  p  ,  densities.  The  approach  of  Reddy  and  Rubin  [12]  where  the 

TD*“  1 

pressure  Is  specified  both  at  inflow  and  at  outflow  is  inconsistent  unless 
one  happens  to  know  the  right  pressures  before  the  computation.  The 
inconsistency  and  consequent  error  can  be  easily  demonstrated  by  one- 
dimensional  examples. 

6.  RESULTS 

In  order  to  check  the  MG  algorithm,  we  choose  the  following  analytic¬ 
al  solution.  It  satisfies  the  continuity  equation  but  gives  rise  to 
source  terms  in  the  momentum  equations: 

U  =  A  +  (x+y)m;  V  =  -(x+y)m;  P  =  - (E1+E2) (x+y ) m  (6.1) 

where  a  and  b  from  equations  (3)  are  defined  by: 

a  =  El  +  F(x+y)n;  b  =  E2  -  F(x+y)n  (6.2) 

and  El  =  1;  E2  =  .  2 ;  F  ■=  .  2 ;  A  =  S ;  Re  =  1000;  m  =  4;  n  =  2.  The 
coarsest  grid  consists  of  4*4  intervals. 

Figure  5  compares  the  MG  convergence  history  of  different  relaxation 
schemes.  In  the  MG  solutions  three  levels  were  involved  (M=3) .  The 
horizontal  coordinate  gives  the  number  of  Work  Units  (WU) ,  where  each 
work  unit  is  equivalent  to  one  global  iteration  on  the  finest  grid.  The 
vertical  coordinate  gives  the  logarithm  of  the  dynamic  residual  c.  The 
dots  show  the  solution  of  the  equivalent  Poisson  equation  (with  the  same 
solution  for  the  pressure  but  with  Dirichlet  condition  over  all  the 
boundaries) .  The  linearized  PNS  equations  were  solved  with  and  without 
the  streamwise  pressure  gradient  correction  of  [3] .  The  corresponding 
(17  *  17  points)  single  grid  convergence  history  is  plotted  for  compar¬ 
ison  (for  the  case  of  Cl  ■  1).  The  corrected  discrete  equations  and 
the  Poisson  equation  exhibit  very  similar  convergence  whereas  the  conver¬ 
gence  of  the  unmodified  equations  is  much  worse.  Upon  increasing  the 
number  of  grids  in  the  unmodified  equations,  the  convergence  deteriorates. 

The  Reynolds  number  independence  of  the  scheme  is  demonstrated  in 
Figure  6,  where  the  convergence  history  is  presented  for  Reynolds  numbers 
1,  10^  and  infinity. 

In  order  to  check  the  non-linear  version  of  the  code,  several  test 


cases  were  run;  the  i nconpr e r r. i bit  flow  over  a  flat  plate,  the*  flow 
along  an  axisymmetric  cylinder,  entrance  flow  between  two  fla*  plates, 
and  the  flow  behind  the  trailing  edge  of  a  flat  plate.  In  all  cases 
good  agreement  was  obtained  with  known  solutions.  The  details  will  be 
presented  elsewhere.  Here  we  show  (Figure  7)  the  convergence  history 
for  a  flow  over  a  flat  plate  with  Uniterm  upstream  profile  and  Neumann 
condition  for  the  pressure  at  the  exit.  K;.ile  the  number  cf  levels  is 
varied, the  finest  grid  remains  the  same  and  consists  of  65  >  65  points. 
In  Figure  8  there  is  a  comparison  between  the  present  results  for  the 
flow  near  the  trailing  edge  of  a  flat  plate  and  the  results  of  refer¬ 
ence  [11].  The  skin  friction  coefficient  CF  is  shown  for  z  <  1  while 
the  center  line  velocity  UC  is  shown  for  z  >  1.  The  trailing  edge  is 
at  z  =  1 . 
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Figure  5:  Convergence  history  for  different  relaxation  schemes  (M  =  3). 
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